Species clustering in competitive Lotka-Volterra models 
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We study the properties of Lotka-Volterra competitive models in which the intensity of the inter- 
action among species depends on their position along an abstract niche space through a competition 
kernel. We show analytically and numerically that the properties of these models change dramat- 
ically when the Fourier transform of this kernel is not positive definite, due to a pattern forming 
instability. We estimate properties of the species distributions, such as the steady number of species 
and their spacings, for different types of kernels. 
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It is widely believed that competition among species 
greatly influences global features of ecosystems. One of 
the most relevant is the fact that ecosystems can host a 
limited number of species. The common explanation is 
the so-called limiting similarity [H and involves represent- 
ing species as points in an abstract niche space, whose 
coordinates quantify the phenotypic traits of a species 
which are relevant for the consumption of resources, usu- 
ally the typical size of individuals, but also preferred prey, 
optimal temperature and so on. On general grounds, one 
expects that a species experiences a stronger competition 
with the closer species in this space. As a consequence, 
a species can survive if it is able to maintain its distance 
with the others above a minimum value which depends 
on the intensity of competition. On the contrary, when 
the distance between two species becomes too small, the 
unavoidable difference in how efflciently the two species 
feed on the resources will make one of the two outcompete 
the other, that will go extinct: this is the essence of the 
competitive exclusion principle [3| , and it can be thought 
as the theoretical grounding of the concept of ecological 
niche. Thus, one expects a stable ecosystem to display a 
flnitc number of species, more or less equidistant in the 
niche space. The finiteness of the number of species has 
been observed in several competition models 3] and, re- 
cently, it has been rigorously demonstrated for a general 
class of them . 

Deviations from the above scenario have aroused re- 
newed interest recently, when it was observed numeri- 
cally [5] that the equilibrium state of rather standard 
models is not always characterized by a homogeneous 
distribution of species in niche space. Instead, clumpy 
distributions, with clusters of many species separated by 
unoccupied regions, are observed. Evidences of a similar 
phenomenon have been observed recently in an evolution- 
ary model [il, suggesting that a theoretical explanation 
of these patterns could bring new insights in the study of 
speciation mechanisms J]. 

In this Letter, we study the Lotka-Volterra (LV) com- 
petitive model as the prototype of competitive systems. 
While the statistical properties of the antisymmetric LV 



model have been addressed recently f§|, not much has 
been demonstrated about the statistics of the competi- 
tive case. Through this paper we will use the language 
of ecological species competition, but we stress that the 
LV set of equations appears in contexts as diverse as 
multimode dynamics in optical systems technology 
substitution [1C| |. or mode interaction in crystallization 
fronts [m , and it is the natural starting point when mod- 
elling competitive systems. Our main result is that the 
macroscopic clustering of species is related to a pattern 
forming transition that separates two different qualita- 
tive behaviors. This is the same phenomenology found 
in birth-death particle systems with interaction at a dis- 
tance, in which individuals typically aggregate forming 
clusters which arrange in an ordered pattern [12], with 
the physical space playing the role of niche space. The 
feature which is relevant for this transition is not the in- 
tensity of the competition, but the functional form of the 
competition kernel. We estimate analytically the num- 
ber of species in cases at both sides of the transition, and 
also discuss the role of species heterogeneity. 

Mathematically, one has competition when the growth 
of a species affects negatively the growth rate of other 
species. One of the simplest implementations of that is 
the LV competitive model: 
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i = 1 ... TV. 
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N is the number of species and Ui denotes the population 
of species i. Each species is characterized by a growth 
rate r and a competition parameter (we take into ac- 
count differences among species only in the latter param- 
eter) . A species is also characterized by a position Xi in a 
niche space that we assume, for simplicity, to be the seg- 
ment [0, L] with periodic boundary conditions. General- 
ization to a multidimensional niche space is straightfor- 
ward, and we expect the unrealistic boundary conditions 
assumed here to be irrelevant except close to the interval 
endpoints. The competition kernel g{x) is a non-negative 
and non-increasing function. Note that the sum in Eq.|T]) 
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contains the self-interaction term g{0)ni. 

To fully specify the dynamics, we should state how the 
Xi are assigned to species and eventually changed. We 
consider a slow immigration rate / at which new species, 
characterized by a random phenotype x € [0,L], are in- 
troduced in the system with a small random population 
6n. This choice is appropriate to model a situation like 
an ecological community on an island [l^ . In addition, 
we consider extinct, and remove from the system, species 
whose population goes below a given threshold ht- When 
is very large compared to the timescales of popula- 
tion dynamics, the system has time to relax to a quasis- 
teady state after each immigration event. Our interest 
here is in the characteristics of these states, in which 
immigration plays almost no role. An efhcient way to 
obtain them is by running simulations of Eq. ([!]) with 
a small amount of immigration, which is then "switched 
off" after some time. By a choice of the time and the 
population units, we can set r — 1 and ht — 1 (thus the 
parameters are really aiUT/r). We also measure all 
distances in niche space in units of L, so that L = 1. 
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FIG. 1: Steady states of system with ai = a = 0.1, 
obtained by evolving a random configuration, initially of 200 
species, for a time t — 5 x 10^. The immigration rate was 
initially / — 0.004, and switched off after half simulation. 
Top panels; Competition kernel gi{x) = exp{~x/ R) (left), 
and g4{x) — exp[—{x/R)'^] (right), with R — 0.1. In the 
bottom panels, we added a Kronecker delta S^o to the kernels 
above. In the last panel, the dotted line is a steady solution 
of ([2]), arbitrarily scaled in the vertical to fit in the same plot. 

In Fig. [1] we show numerical results for the dynam- 
ics just defined, for Ui — a for all i. In the top row we 
compare the distribution of species with kernels gi{x) — 
exp(— cc/i?) and g4{x) ~ cxp[— (a;/i?)^]. In the exponen- 
tial case (left) species occupy the full niche space. Al- 
though they are not perfectly equispaced and there are 
differences in the population sizes, there is a clear aver- 
age interspecies distance, which corresponds to 1/N. In 
the quartic-exponential case (right), a much more regular 



pattern emerges, with different species perfectly equidis- 
tant (and all with the same population). Since growth 
limitation is known to affect species distribution [5|, we 
plot in the bottom row results for kernels gi (x) -|- S^o and 
54 (x) -|- 6x0, be. the same kernels but with an enhanced 
value of the self-competition coefficient g(0). Here the 
difference is even more striking: the exponential case is 
similar to the previous one, but the quartic case shows 
clear clusters of species separated by empty regions. 

To understand the origin of the periodic patterns, we 
write a continuum evolution equation for the field (f>{x, t), 
the expected density of individuals in a given point x of 
the niche space as a function of time: 

dt<p{x, t) = (j){x, t){l-a / g{\x - y\)(j){y, t)dy ] + s, 



(2) 

which is a mean field version of Eq. ([1]) for — a. 
In this macroscopic description, we neglect fluctuations 
in the immigration process by using a constant rate 
s — ISn. The stationary homogeneous solutions of 
Eq. (HD are (po = {1 ± VT+Tsa) / {2a) , where a = oAf, 
with J\f — J g{x)dx. Of the two solutions, only the one 
with the plus sign is acceptable since the other leads to 
a negative density (this second solution corresponds to 
the extinct absorbing state when s = 0). We analyze the 
stability of the positive solution by considering a small 
harmonic perturbation </) = 00 + eexp(Ai -I- ikx). Sub- 
stituting into ([2]), the first order in e gives the following 
dispersion relation: 
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where g{k) = J 9{^) exp{—ikx)dx is the Fourier trans- 
form of g{x). When A becomes positive for some values 
of k, the constant solution of© is unstable, signaling 
a pattern forming transition [14| with the characteristic 
length scale of the pattern determined by the value of 
k at which A(fc) is maximum. For Eq. ([3]), in the limit 
s ^ 0, it is a sufficient and necessary condition for insta- 
bility that the Fourier transform of the kernel, ^, takes 
negative values (notice that (j)Qd > 1, with (/jqci ^ 1 in 
the limit s — > 0; for sufficiently large s, the homogeneous 
state regains stability). To exemplify this mechanism, 
we consider the family of kernels ga{x) = exp[— (x/i?)'^], 
being a > and R the typical competition range. It is 
known [151] that this family of functions has non-negative 
Fourier transform for < a < 2. Interestingly enough, 
the Gaussian kernel, which is the commonly adopted one 
[J, 01, corresponds to the marginal case. This may imply 
that some results previously obtained for this case could 
be non robust and largely affected by the way immigra- 
tion is introduced, the presence or absence of diffusion 
processes in niche space, etc. 

We quantify the pattern-forming transition in terms 
of the structure function, S{k) — \ rij e'xjp{ikxj)\'^ , of 
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the stationary distribution of species obtained from the 
simulations. The position and height of its maximum 
identify periodic structures. In Fig. Owe plot (left panel) 
the maximum height of S* as a function of the exponent a 
of the kernel. The sharp increase of max S for a > 2 in- 
dicates the formation of periodic structures in this range. 
This is confirmed by the right panel plot, where we show 
the position km of the peak of S, together with the value 
fci at which the linear growth, expression Q, has a max- 
imum. Note that the location of this maximum is inde- 
pendent of the parameters a and s, being only dependent 
on the parameters in ga{x) {R and cr; the dependence 
on R disappears when considering k^R). The striking 
agreement between km and k^ for a > 2 confirms that 
the linear pattern forming instability of the homogeneous 
distribution is the mechanism responsible for the peri- 
odic species arrangement observed in that range. Except 
when (T 2, the value of k^R is in the range 4.5-5.0, so 
that the pattern periodicity would be d « 27r/fci « aR, 
with a « 1.3 — 1.4, as observed in Fig. [1] (right panels). 

Another difference between a < 2 and a > 2, visible 
in Fig. [U is the existence in the later case of exclusion 
zones around established species, in which immigrants 
have not been able to settle. We can understand the 
presence of these regions also from the density equation 
([2]), for s = 0, by noticing that its steady stable solu- 
tions (j)st{x) necessarily have regions with (pstix) = in 
the pattern forming case. This can be seen from the 
steady state condition J dyg{\x — y\)(f>st{y) = 1/a, which 
is valid at all niche locations x at which (pstix) ^ 0. If 
these locations cover in fact the full niche space [0, 1], we 
can solve the steady state condition by Fourier transform 
and find that the only solution (for nonconstant g{x)) is 
the homogeneous one 4>st{x) — {aj\f)~^. Since we know 
that this is linearly unstable when the Fourier transform 
of ^(a;) is not positive definite, we conclude that steady 
stable solutions of ([2|) in the pattern forming case must 
have regions of zero density, which we identify with the 
exclusion zones. Given the absorbing character of the 
(j) — state, many steady solutions exist, differing in 
the amount and location of the (f>st = segments, but 
the most relevant are the ones attained when s — > 0+. 
Figure [1] (bottom right) shows one of these solutions, nu- 
merically obtained (for a kernel g4{x)+S{x)). The steady 
solution corresponding to the g4{x) kernel of the top right 
panel is zero everywhere except at a set of periodically 
spaced delta functions. In both cases the discrete species 
distribution is well represented by the solutions of ([2]). 

When g{k) remains positive, as for go{x) with cr < 2, 
A(fc) remains negative, and there are no patterns nor ex- 
clusion zones surviving in steady solutions of the density 
equation for s — > 0+. Thus, the characteristic distance 
between species, observed in Figs. [Hand [2] to be quahta- 
tively different from the case a > 2, should be determined 
by a different mechanism. We explore it for the exponen- 
tial kernel, gi{x) = exp(— x/i?), because it allows some 
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FIG. 2: Left panel: maximum of {S{k)), the structure func- 
tion averaged over 1000 realizations of stationary distributions 
of species (obtained after a time t = 10"", without immigra- 
tion during the last half of it) as a function of a, for R = 0.1, 
a = 0.1. Right panel: position of the peak km vs a (circles), 
together with the linearly fastest growing mode A;l (line), from 
([3]). For (T > 2, the difference between km and /cl is always 
smaller than the finite-size discretization of the values of km- 



analytical estimates. Fig. [3] shows the number of species 
at equilibrium for — a, and also in the heterogeneous 
situation in which the Ui 's are independent random vari- 
ables uniformly distributed between 0.95a and 1.05a, be- 
ing a an average value. 




FIG. 3: Number of species as a function of competition co- 
efficient a (left panel, i? = 0.1) and the interaction distance 
R (right panel, a — 0.1). Symbols joined by dashed and 
solid lines are for the cases of a'^s heterogeneity (for which a 
is plotted instead of a), and non- heterogeneity, respectively. 
The upper dot-dashed lines are from the approximation (|4]). 



During the system evolution (in the non-heterogeneous 
case) we observe that, when the species are far from 
the extinction threshold, it is always possible for a new 
species to successfully settle between two of them. But 
as a consequence the populations of these neighboring 
species get reduced. This brings the species populations 
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closer to ht as the number of species increases, and even- 
tually no new species will be admitted. Thus, in this 
a < 2 case, the mechanism fixing a maximum number of 
species, and thus a characteristic mean distance d among 
them, is the presence of the extinction threshold ut- 

Fig. [3] shows that the number of species N (in the 
steady state obtained after switching off immigration) 
grows linearly with 1/ R and with 1/a in the case of equal 
species, while heterogeneity slows down the increase with 
1/R and almost stops it with 1/a. Notice that decreasing 
a is the same as decreasing the threshold value ut due 
to our rescaling of the equations. We can explain these 
dependences by considering an ideal steady state made 
of equidistant species, at distance d, and having the same 
population n* . The equilibrium condition for the system 
of equations ([1]) in the exponential kernel case becomes 
tanh((i/2i?) — an*, which gives d « 2aRn* in the limit 
{d/2R) ^ 1. Recalling that N — 1/d, each population 
n* decreases as the number of species increases during 
immigration. The limit, setting the steady state, will be 
the situation in which n* = ut = 1, for which no new 
immigrant can be accepted. Thus we estimate the equi- 
librium number of species in this case as 

N « {2aRy^ . (4) 

This is only a rough approximation, since species are 
not equidistant nor equipopulated in the true equilib- 
rium, but provides an explanation for the observed linear 
scaling of N with 1/a and 1/R. Fig. [3]shows that it gives 
an upper estimation for the number of species in the less 
ordered distributions actually found. 

The case with heterogeneity of Fig. [3] shows a clearly 
different mechanism: the number of species does not 
change with a and consequently with ht. Thus, this 
scenario is qualitatively similar to the pattern forming 
case: there is a distance in the niche space, not related 
to the threshold value, of the order of the interaction 
range. Two species cannot survive due to heterogeneity 
if they are closer than this distance, independently on 
the mean competition strength. 

To clarify this mechanism, we consider the role of het- 
erogeneity in the long-range case of a constant kernel, 
g{x) = 1 for all x. This may be interpreted as a case in 
which the kernel decaying distance goes to infinity. Sum- 
ming all the equations in ([T]) we obtain an equation for 
the total population Ntot — J2j "-j • 

Ntot = Ntot{l - {a)Ntot) (5) 

where (a) = (/^j a,jnj)/Ntot- After a short time, the 
equilibrium value Ntot = i^-)^^ would be attained and 
we can plug this value back into Eqs. ([T]) to obtain 
"hi = ni(l — ai/{a)), valid at longer times. In the case of 
equal species, one has ai = a = (a) and all possible states 
with a rij = 1 are allowed. In the heterogeneous case, 
species having < (a) will grow while the others will 



decrease their population and finally go extinct. Mean- 
while, it is easy to realize that (a) will increase, sending 
more and more species below the extinction threshold. 
The final result, valid for any initial distribution of the 
fli's, is that just one species will survive, as confirmed by 
simulations (not shown). 

To conclude, we studied analytically and numerically 
the collective behavior of competitive Lotka-Volterra sys- 
tems. Our main message is that the form of the competi- 
tion kernel changes drastically the equilibrium distribu- 
tion of species. Species clustering with periodic spacings 
of the order of the interaction range can occur at one side 
of a pattern forming transition, whereas smaller spac- 
ings, depending on the interaction strength a, occur at 
the other. Surprisingly, the Gaussian kernel, the one usu- 
ally considered in the literature, corresponds to a fron- 
tier case. Diversity has been shown to alter qualitatively 
the competition outcome. Generalization to the case of a 
multidimensional niche space is straightforward and does 
not change qualitatively our results. Diffusion in niche 
space, modelling mutations 0], can be introduced and 
has a stabilizing effect somehow similar to that of the 
immigration rate s. 
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